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Abstract 

Background: Ion channels play a critical role in a wide variety of biological processes, including the development of 
human cancer. However, the overall impact of ion channels on tumorigenicity in breast cancer remains controversial. 

Methods: We conduct microarray meta-analysis on 280 ion channel genes. We identify candidate ion channels that 
are implicated in breast cancer based on gene expression profiling. We test the relationship between the expression of 
ion channel genes and p53 mutation status, ER status, and histological tumor grade in the discovery cohort. A 
molecular signature consisting of ion channel genes (IC30) is identified by Spearman's rank correlation test conducted 
between tumor grade and gene expression. A risk scoring system is developed based on IC30. We test the prognostic 
power of IC30 in the discovery and seven validation cohorts by both Cox proportional hazard regression and log-rank 
test. 

Results: 22, 24, and 30 ion channel genes are found to be differentially expressed with a change in p53 mutation 
status, ER status, and tumor histological grade in the discovery cohort. We assign the 30 tumor grade associated ion 
channel genes as the IC30 gene signature. We find that IC30 risk score predicts clinical outcome (P < 0.05) in the 
discovery cohort and 6 out of 7 validation cohorts. Multivariate and univariate tests conducted in two validation 
cohorts indicate that IC30 is a robust prognostic biomarker, which is independent of standard clinical and pathological 
prognostic factors including patient age, lymph node status, tumor size, tumor grade, estrogen and progesterone 
receptor status, and p53 mutation status. 

Conclusions: We identified a molecular gene signature IC30, which represents a promising diagnostic and prognostic 
biomarker in breast cancer. Our results indicate that information regarding the expression of ion channels in tumor 
pathology could provide new targets for therapy in human cancers. 
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Background immune response, cell volume regulation, and cell prolif- 
Ion channels are membrane proteins expressed in vari- eration [1-7]. Because of the involvement of ion channels 
ous tissues that allow the passage of ions across bio- in diverse biological functions, defects in the expression 
logical membranes. Ion transport is a key component in and functional activity of ion channels can cause disease 
a wide variety of biological processes including electrical in many tissues [8]. The number of human diseases re- 
impulse generation and conduction along nerves, fluid lated with ion channel malfunction has grown rapidly 
balancing within cells and across cell membranes, and over the past few years [4,9,10]. In particular, there is in- 
signal transduction within and among cells. In addition, creasing evidence that ion channels, including both 
ion channels are known to play critical roles in gene voltage-gated and ligand-gated channels, are involved in 
expression, hormone secretion, muscle contraction, the progression and pathology of diversified human can- 
cers [6,7,11-17]. For example, voltage-gated potassium 
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ation through the modulation of membrane potential in 
breast, colon, and prostate cancers [12,14,15]. Transient 
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receptor potential (TRP) channels are involved in vascular 
permeability and angiogenesis and have been implicated in 
tumor growth and metastasis [18,19]. Several ligand-gated 
channels, such as nicotinic acetylcholine receptors, affect 
neoplastic progression by regulating tumor cell prolifera- 
tion, apoptosis, and angiogenesis [20-23]. More impor- 
tantly, the expression level of ion channels is potentially 
able to serve as a prognostic index in human cancers 
for clinical purposes. For instance, TRPMl, a TRP cat- 
ion channel, is an indicator of melanoma aggressiveness 
[24] and expression of the Ca^"^-selective cation channel 
TRPV6 is a prognostic marker for tumor progression in hu- 
man prostate cancer [25]. In addition, the long TRP channel 
TRPM8 might serve as a prognostic marker for androgen- 
unresponsive and metastatic prostate cancer [26] and 
the expression of SCN9A, a voltage-gated sodium (Na"^) 
channel, is also useful for prognostic purposes in pros- 
tate cancer [27]. 

Breast cancer is the most common invasive cancer in 
women worldwide. It is also the principal cause of death 
from cancer among women globally [28]. A large num- 
ber of breast cancer studies sought to understand the 
molecular mechanisms of cancer origin, progression, 
and invasion that lead to metastasis, and many of these 
studies have underlined the involvement of ion channels 
in breast cancer. For example, Kca channels contribute 
to breast tumor migration and progression [29,30] and 
TRP channels are strongly correlated with breast tumor 
cell proliferation [31]. In addition, the upregulation of 
several voltage-gated Na^ (NaV) channels is associated 
with metastatic process in breast cancer [32]; however, 
most of these studies focus on only one ion channel or 
one type of ion channel. So far, there is no clear picture 
on the overall expression profiling of different ion chan- 
nel genes in breast cancer. High-throughput "omic" 
technologies make it possible to scan all ion channel 
genes rather than focusing on a single gene or gene fam- 
ily [33]. In this study, we looked to identif)^ molecular 
signatures consisting of multiple genes from different 
ion channel families that are implicated in the pathology 
of human breast cancer. Firstly, we investigated the asso- 
ciation of ion channel genes with p53 mutation status in 
breast tumors. The tumor suppressor p53 is known to 
play a critical role in regulating the cell cycle and is thus 
involved in preventing cancer. Mutations in p53 are 
strongly associated with poor clinical outcome in breast 
cancer patients [15,23]. Comparison between the p53 
mutant and wild-type groups showed that ion channel 
genes are associated with more aggressive and thera- 
peutically refractory tumors [15]. Secondly, we identified 
the ion channel genes that were differentially expressed 
between estrogen receptor (ER) positive and negative 
breast cancer patients. About 75% of all breast cancers 
are ER positive, which grow in response to the hormone 



estrogen. ER is a powerful prognostic marker and an ef- 
ficient target for the treatment of hormone-dependent 
breast cancer [26]. Identification of the ER-related ion 
channels helps us understand the role of ER in the de- 
velopment and progression of breast cancer. Thirdly, we 
investigated the relationship between ion channel gene 
expression and histological tumor grade in breast cancer. 
We identified a molecular signature consisting of 30 ion 
channel genes (IC30), which significantly correlated with 
tumor grade. We demonstrate that IC30 is a robust 
prognostic biomarker to predict clinical outcome in 
breast cancer, and is independent of standard clinical 
and pathological prognostic factors including patient 
age, lymph node status, tumor size, tumor grade, ER sta- 
tus, and progesterone receptor (PR) status. The perform- 
ance of IC30 was also validated in several independent 
cohorts from different parts of the world (Table 1). 

Results 

Differentially expressed ion channel genes between p53 
mutant and wild-type tumors 

280 ion channel genes were collected for this study 
(Additional file 1: Table SI). We aim to identify candi- 
date ion channels that are implicated in breast cancer 
based on gene expression profiling. We first explored 
the difference in ion channel gene expression between 
p53 mutant and wild-type breast tumors in the discovery 
SIN cohort. There were 58 samples with p53 mutations 
resulting in protein-level changes and 193 samples with 
a wild- type genotype [15]. In total, 22 ion channel genes 
were identified as differentially expressed between 
the two groups. Five ion channel genes were upregulated 
in p53 mutant tumors, including KCNE3, KCNN4, 
and MC0LN2; while 17 ion channel genes were 
downregulated, including ANOl, KCNMAl, and TPCNl 
(Table 2 and Figure 1). Among these differentially 
expressed genes, all the Ca^^ channel {CACNAID, 
CACNA2D1, and CACNA2D2) and Na"" (.SCA^Z4 and 
SCNNIA) channel genes were downregulated in mutant 
tumors. In contrast, the expression pattern of chan- 
nel and chloride (CI ) channels was more heterogeneous. 
Genome-scale inspection indicated a significant enrich- 
ment of ion channel genes among the genes regulated 
by p53 mutation status {P = 0.027 by Fisher s exact test). 

To test the reliability of the above results in another 
cohort, we accessed a publicly available microarray 
dataset on breast cancer (FRA) where p53 mutation sta- 
tus was known. Unsupervised hierarchical cluster ana- 
lysis demonstrated a very similar expression pattern of 
the 22 differentially expressed ion channel genes be- 
tween the SIN and FRA cohorts (Figure 2). We evalu- 
ated the statistical significance in hierarchical cluster 
analysis by approximately unbiased P- value {AU) (see 
Methods for details). In the SIN cohort, the hierarchical 
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Table 1 Gene expression datasets of breast cancer 



Organization of data source 


Abbreviation 


GEO accession 


Reference 


(^pnnmp In^^titiitp nf '^in^;^^^^p '^in^;^^^^p 
vjci lui 1 ic II loLiLULC wi jii lyapwic, jii iyci|-j>>jic 


SIN 


GSE3494 




Institut Paoli-Calmettes Marseille, France 


FRA 


GSE21653 


[16] 


Siemens Medical Solutions Diagnostics GmbH, Germany 


GER 


GSE11121 


[34] 


Netherlands Cancer Institute, Netherlands 


NED 


a 


[21] 


Karolinska Instituted Sweden 


SWE 


GSE1456 


[18] 


Koo Foundation SYS Cancer Center, Taiwan 


7WN 


GSE20685 


[19] 


Nuvera Biosciences Inc, United States 


USA1 


GSE25066 


[20] 


Veridex LLC, United States 


USA2 


GSE2034 


[27] 



^Expression data for NED are publicly available from http://bioinformatics.nki.nl/data.php. 



clusters of upregulated and downregulated genes were 
highly robust {AU= 0.993 for the upregulated cluster 
and AU= 0.963 for the downregulated cluster). Similar 
results were obtained for the FRA cohort {AU = 0.985 
for the upregulated cluster and AU= 0.990 for the 
downregulated cluster). Two-tailed t- tests also indicated 
that 15 out of the 22 genes were significantly (adjusted 
P < 0.05) dysregulated between p53 mutant and wild- 
type tumors in the FRA cohort (Additional file 1: Table 



S2). The direction of differential expression in the SIN 
cohort was reproduced in the FRA cohort (Figure 3A). 

Differentially expressed ion channel genes between ER 
positive and negative patients 

We compared ion channel gene expression between ER 
positive and negative patients in the SIN cohort. A total 
of 213 patients were identified as ER positive while 34 
patients were identified as ER negative. Twenty-four ion 



Table 2 Comparison in gene expression level between p53 mutant and wild-type tumors 



Gene 
symbol 


Gene title 


Fold 
change^ 


Adjusted 
P-value'' 


ANOl 


anoctamin 1, calcium activated chloride channel 


0.66 


0.001 


CACNAID 


calcium channel, voltage-dependent, L type, alpha ID subunit 


0.32 


< 0.001 


CACNA2D1 


calcium channel, voltage-dependent, alpha 2/delta subunit 1 


0.79 


0.002 


CACNA2D2 


calcium channel, voltage-dependent, alpha 2/delta subunit 2 


0.47 


< 0.001 


CLCA2 


chloride channel accessory 2 


2.53 


0.039 


cues 


chloride intracellular channel 5 


0.77 


< 0.001 


CUC6 


chloride intracellular channel 6 


0.28 


< 0.001 


GLRB 


glycine receptor, beta 


0.32 


< 0.001 


KCND3 


potassium voltage-gated channel, Shal-related subfamily, member 3 


0.46 


< 0.001 


KCNE3 


potassium voltage-gated channel, Isk-related family, member 3 


1.34 


< 0.001 


KCNE4 


potassium voltage-gated channel, Isk-related family, member 4 


0.25 


< 0.001 


KCNJ3 


potassium inwardly-rectifying channel, subfamily J, member 3 


0.51 


< 0.001 


KCNKl 


potassium channel, subfamily K, member 1 


1.44 


< 0.001 


KCNK6 


potassium channel, subfamily K, member 6 


0.80 


0.010 


KCNMAl 


potassium large conductance calcium-activated channel, subfamily M, alpha member 1 


0.65 


< 0.001 


KCNN4 


potassium intermediate/small conductance calcium-activated channel, subfamily N, member 4 


1.69 


< 0.001 


MC0LN2 


mucolipin 2 


1.60 


< 0.001 


P2RX4 


purinergic receptor P2X, ligand-gated ion channel, 4 


0.73 


< 0.001 


SCN7A 


sodium channel, voltage-gated, type VII, alpha subunit 


0.29 


< 0.001 


SCNNIA 


sodium channel, non-voltage-gated 1 alpha subunit 


0.71 


< 0.001 


TPCNl 


two pore segment channel 1 


0.77 


< 0.001 


TRPCl 


transient receptor potential cation channel, subfamily C, member 1 


0.72 


0.005 



^Fold change is calculated by dividing the mean expression of p53 mutant tumor by the mean expression of p53 wild-type tumor. 
■^P-value is calculated by two-tailed f-test and adjusted by Benjamini & Hochberg correction. 
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Figure 1 Boxplot of expression of the genes differentially expressed between p53 mutant and wild-type tumors. Twenty-two ion 
channel genes were found to be differentially expressed between the two groups. Light-grey represents wild-type group while dark-grey 
represents mutant group. Y-axis: log2-transformed expression values. 



channel genes were identified as differentially expressed 
between the two groups; 16 genes were upregulated in 
ER positive patients while 8 genes were downregulated 
(Table 3 and Figure 4). Among these differen- 
tially expressed genes, all the Ca^"^ channels (except 
CACNAIA) and Na"^ channel genes were upregulated in 
the ER positive group, whereas the expression pattern of 
K"^ channel and CI' channels was more heterogeneous. 
Nineteen out of the 24 differentially expressed genes 
overlapped with the genes differentially expressed be- 
tween p53 mutant and wild-type tumors (Figures 1 and 
4), which is statistically significant (P < 0.001 by cumula- 
tive hypergeometric distribution function). Among these 
common genes, all the downregulated genes in ER 
positive patients were upregulated in the p53 mutant 
group and vice versa, which is consistent with previous 
findings that patients with negative ER status share com- 
mon pathology with patients harboring mutant p53 
[25,35]. We also found a significant enrichment of ion 
channel genes among the genes dysregulated by ER sta- 
tus (P = 0.013 by Fisher s exact test). 



We next checked the expression profiling the 24 ER 
status related ion channel genes in three independent 
cohorts (FRA, USAl, and USA2) where ER status was 
known. The heatmaps generated by unsupervised hier- 
archical cluster analysis demonstrated an analogous ex- 
pression profiling for the 24 ion channel genes among 
the discovery and validation cohorts (Additional file 2: 
Figures SI, S2, S3 and S4). A side-by-side comparison 
between Table 3 and Additional file 1: Table S3 revealed 
that the 24 genes were significantly (adjusted P < 0.05) 
dysregulated by ER status in at least one out of the three 
validation cohorts. The direction of differential expres- 
sion in the SIN cohort was consistent with that in the 
FRA, USAl, and USA2 cohorts (Figure 3B). 

Correlation between tumor grade and expression of ion 
channel genes 

To determine the relationship between tumor progres- 
sion and ion channel gene expression, we linked gene 
expression level with histological tumor grade in the SIN 
cohort, using the Spearman's rank correlation test. The 
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Figure 2 Heatmaps of expression of the genes differentially expressed between p53 mutant and wild-type tumors. The differentially 
expressed ion channel genes were derived from the discovery cohort (SIN) and verified in the FRA cohort. Each row in the heatmaps was 
labelled with the corresponding gene symbol. The columns labelled with "-" denote p53 mutant tumors. Red represents relatively increased gene 
expression while blue represents down-regulation. The horizontal dotted line separates the down- and up- regulated gene clusters. 



expression of 30 ion channel genes was found to be sig- 
nificantly (adjusted P < 0.05) correlated with tumor grade 
(Table 4). Eleven out of the 30 genes showed a positive 
correlation between expression and tumor grade while 
the other 19 genes showed a negative pattern (Figure 5). 
Given the fact that tumor grade reflects the differenti- 
ation of breast cancer cells, we identified 11 upregulated 
and 19 downregulated ion channel genes in more 
aggressive breast tumors. Among these 30 genes, 19 
genes were also differentially expressed between the p53 
mutant and wild-type tumors. The number of the over- 
lapping dysregulated genes was statistically significant 
(P< 0.001 by cumulative hypergeometric distribution 
function). Positive correlation between expression and 
tumor grade corresponds to upregulation in p53 mutant 
tumors and vice versa, which confirms the well- 
established findings that p53 mutations link to higher- 
grade breast cancer and potentially poorer clinical 
outcomes [32,33,36,37]. Similar to the results described 



above for p53/ER status, ion channel genes were also 
significantly enriched among all the genes that were re- 
gulated by tumor grade {P = 0.003 by Fisher s exact test). 

Because of the availability of tumor grade information 
in the FRA, GER, and USAl cohorts, we tested the ex- 
pression pattern of the above 30 genes in these 3 inde- 
pendent datasets. We observed a significant correlation 
between the expression and tumor grade in at least 
one validation cohort for each gene, except for TPCN2 
(Additional file 1: Table S4). The correlation coefficients 
for each gene were largely concordant across the discovery 
and validation cohorts (Figure 3C). 

Prognostic molecular signature composed of ion 
channel genes 

We hypothesized that the 30 ion channel genes associ- 
ated with tumor histological grade would be predictive 
of tumor outcome in breast cancer patients. We 
designated these ion channel genes as the IC30 gene 
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Figure 3 Validation of ion ciiannel gene expression profiling in breast cancer. (A) Heatmap of log2-transfornned fold change of gene 
expression level between p53 mutant and wild-type tumors (mutant/wild-type). In the SIN cohort, 5 genes were up-regulated (log2-transformed 
fold change > 0) while 17 genes were down-regulated (log2-transformed fold change < 0) in mutant tumors. The vertical dotted line separates 
the down- and up- regulated gene clusters {AU= 1 for both clusters). (B) Heatmap of log2-transformed fold change of gene expression level 
between ER positive and negative groups (ER positive/negative). In the SIN cohort, 16 genes were up-regulated (log2-transformed fold change > 
0) while 8 genes were down-regulated (log2-transformed fold change < 0) in ER positive patients. The white cells in the heatmap mean the gene 
expression data are unavailable in the corresponding datasets. The vertical dotted line separates the down- and up- regulated gene clusters {AU = 
0.995 for both clusters). (C) Heatmap of Spearman's rank correlation coefficient. We found 30 ion channel genes with significant Spearman's rank 
correlation between gene expression and histological grade in SIN cohort. Correlation coefficient > 0 means that gene is upregulated in the patients 
with higher grade while negative correlation coefficient indicates down-regulation in the patients with higher grade. The white cells in the heatmap 
mean the gene expression data are unavailable in the corresponding datasets. The vertical dotted line separates the down- and up- regulated gene 
clusters {AU= 1 for both clusters). 



signature. We developed a risk scoring system that com- 
bined gene expression information in the IC30 with the 
Spearman's ranl< correlation coefficients (p) listed in 
Table 4. IC30 positive patients were defined as those 
having a risk score greater than the group median score 
and the other patients were assigned as IC30 negative. 



We tested the ability of the risk score to stratify patients 
into prognostic groups in the SIN cohort and the seven 
validation cohorts (FRA, GER, NED, SWE, TWN, USAl, 
and USA2). Kaplan-Meier survival analyses were used to 
compare the IC30 positive and negative groups and de- 
termined statistical significance by log-rank tests. The 
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Table 3 Comparison in gene expression level between ER positive and negative tumors 





VJCI IC IIIIC 


Fold 


AHii ictpH 


symbol 




change^ 


P-value'^ 


ANOl 


anoctamin 1, calcium activated cliloride cliannel 


2.09 


< 0.001 


CACNAIA 


calcium channel, voltage-dependent, P/Q type, alpha 1 A subunit 


0.65 


0.043 


CACNAID 


calcium channel, voltage-dependent, L type, alpha ID subunit 


5.59 


< 0.001 


CACNA2D1 


calcium channel, voltage-dependent, alpha 2/delta subunit 1 


1.30 


0.033 


CACNA2D2 


calcium channel, voltage-dependent, alpha 2/delta subunit 2 


2.61 


< 0.001 


CLCA2 


chloride channel accessory 2 


0.46 


0.017 


CLIC4 


chloride intracellular channel 4 


0.75 


0.008 


CLIC6 


chloride intracellular channel 6 


5.86 


< 0.001 


GABRP 


gamma-aminobutyric acid (GABA) A receptor, pi 


0.20 


0.020 


GLRB 


glycine receptor, beta 


3.42 


< 0.001 


KCNAB2 


potassium voltage-gated channel, shaker-related subfamily, beta member 2 


0.72 


0.013 


KCND3 


potassium voltage-gated channel, Shal-related subfamily, member 3 


2.36 


< 0.001 


KCNE3 


potassium voltage-gated channel, Isk-related family, member 3 


0.75 


0.003 


KCNE4 


potassium voltage-gated channel, Isk-related family, member 4 


12.64 


< 0.001 


KCNJ3 


potassium inwardly-rectifying channel, subfamily J, member 3 


6.22 


< 0.001 


KCNK6 


potassium channel, subfamily K, member 6 


1.82 


< 0.001 


KCNMA 1 


potassium large conductance calcium-activated channel, subfamily M, alpha member 1 


1 .45 


0.035 


KCNN4 


potassium intermediate/small conductance calcium-activated channel, subfamily N, member 4 


0.42 


< 0.001 


KCNS3 


potassium voltage-gated channel, delayed-rectifier, subfamily S, member 3 


1.58 


0.028 


MC0LN2 


mucolipin 2 


0.46 


0.002 


P2RX4 


purinergic receptor P2X, ligand-gated ion channel, 4 


1.82 


< 0.001 


SCN7A 


sodium channel, voltage-gated, type VII, alpha subunit 


1.71 


0.009 


SCNNIA 


sodium channel, non-voltage-gated 1 alpha subunit 


2.00 


< 0.001 


TPCNl 


two pore segment channel 1 


1.38 


< 0.001 



^Fold change is calculated by dividing the mean expression of ER positive tumor by the mean expression of ER negative tumor. 
■^P-value Is calculated by two-tailed f-test and adjusted by Benjamlnl & Hochberg correction. 



IC30 signature was able to identify patients with poor 
breast cancer survival in all the cohorts {P < 0.01) except 
for the GER cohort (a marginal P = 0.066) (Figure 6). 
This association between IC30 status and survival was 
also confirmed by univariate Cox proportional hazard 
regression of survival. IC30 positive patients had a 1.98- 
fold increased risk for death in the SIN cohort, 1.99-fold 
in the FRA cohort, 1.73-fold in the GER cohort, 1.81- 
fold in the NED cohort, 4.33-fold in the SWE cohort, 
1.82-fold in the TWN cohort, 3.11-fold in the USAl co- 
hort, and 1.71-fold in USA2 cohort (Table 5). 

Independence of IC30 from other clinlcopathological 
factors 

We investigated the performance of the IC30 signature 
in comparison with clinlcopathological variables associ- 
ated with prognosis in breast cancer in the USAl cohort, 
the largest dataset in this study. A multivariate Cox re- 
gression of survival indicated that IC30 status remained 
a significant covariate in relation to the standard 



clinlcopathological factors in breast cancer, including 
patient age, lymph node status, tumor size, tumor grade, 
and ER and PR status (Table 6). Patients were stratified 
according to respective clinlcopathological factors. For 
patients aged < 50 and > 50, the IC30 positive patients 
had a significant 2.37-fold (P = 0.003) and 4.21-fold {P < 
0.001) increased risk for death, respectively. For patients 
with and without lymph node involvement, the IC30 
positive patients had a 2.05-fold (P = 0.157) and 2.72-fold 
{P < 0.001) increased risk for death, respectively. For pa- 
tients with tumor size < T3 and > T3, the IC30 positive 
patients had a significantly increased risk for death of 
3.61-fold (P< 0.001) and 2.78-fold (P< 0.001), respect- 
ively. For patients with lower (1 or 2) and higher (3) 
tumor grade, the IC30 positive patients had a signifi- 
cantly 6.91-fold {P < 0.001) and 1.67-fold {P = 0.044) in- 
creased risk for death, respectively. For patients with ER 
negative and positive status, the IC30 positive patients 
had an increased risk for death of 1.30-fold {P = 0.275) 
and 2.94-fold (P = 0.002), respectively. Lastly, IC30 
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Figure 4 Boxplot of expression of the genes differentially expressed between ER positive and negative patients. Twenty-four ion 
channel genes were found to be differentially expressed between the two groups. Light-grey represents ER negative group while dark-grey 
represents ER positive group. Y-axis: log2-transformed expression values. 



positive patients with PR negative and positive status 
had a significantly 1.65-fold (P = 0.030) and 2.35-fold 
(P = 0.021) increased risk for death, respectively. Kaplan- 
Meier survival analysis also demonstrated a significantly 
reduced survival for IC30 positive patients in each sub- 
set grouped by age, lymph node status, and tumor size 
(Figure 7). In addition, univariate Cox regressions of sur- 
vival confirmed that the IC30 signature was a superior 
survival predictor in the USAl cohort, in comparison 
with age, tumor size, and tumor grade (Additional file 1: 
Table S5). 

We also checked the independent predictive power of 
the IC30 signature in the FRA cohort, where informa- 
tion on age, tumor grade, ER status, PR status, and p53 
mutation status is available. A multivariate Cox regres- 
sion of survival indicated that the IC30 signature 
was the only significant survival predictor (P = 0.014) 
(Table 7). Univariate Cox regressions of survival also 
confirmed that the IC30 signature was the most signifi- 
cant prognostic factor in the FRA cohort, in comparison 



with age, tumor grade, ER and PR status, and p53 muta- 
tion status (Additional file 1: Table S6). Taken together, 
these results suggest that IC30 is associated with clinical 
outcome and is an independent prognostic factor. 

Discussion 

Ion channels are implicated in many physiological pro- 
cesses and also play a pivotal role in the development of 
cancers; however, it is currently difficult to assign a spe- 
cific mechanism for each ion channel in the prolifera- 
tion, invasion, and metastasis of tumor cells [17,24]. 
Here, we investigated the pathological role of ion chan- 
nel genes in breast cancer with respect to gene expres- 
sion level. We tested the association of ion channel 
genes with p53 mutation status, ER status, and tumor 
histological grade: 22 ion channel genes were found to 
be dysregulated between p53 mutant and wild-type tu- 
mors, 24 ion channel genes were differentially expressed 
between ER positive and negative patients, and the ex- 
pression level of 30 ion channel genes was significantly 



Ko et al. Molecular Cancer 201 3, 1 2:1 06 
http://www.nnolecular-cancer.conn/content/1 2/1/1 06 



Page 9 of 1 7 



Table 4 Correlation between gene expression and histological tumor grade 



Gene 


Gene title 




Adjusted 


symbol 






P-value*" 


ANOl 


anoctamin 1, calcium activated cliloride cliannel 


-0.23 


< 0.001 


CACNAID 


calcium channel, voltage-dependent, L type, alpha 1 D subunit 


-0.42 


< 0.001 


CACNA2D1 


calcium channel, voltage-dependent, alpha 2/delta subunit 1 


-0.28 


< 0.001 


CACNA2D2 


calcium channel, voltage-dependent, alpha 2/delta subunit 2 


-0.30 


< 0.001 


CLIO 


chloride intracellular channel 1 


0.26 


< 0.001 


CUC4 


chloride intracellular channel 4 


0.16 


0.022 


cues 


chloride intracellular channel 5 


-0.22 


0.001 


CUC6 


chloride intracellular channel 6 


-0.33 


< 0.001 


GLRB 


glycine receptor, beta 


-0.35 


< 0.001 


KCNAB2 


potassium voltage-gated channel, shaker-related subfamily, beta member 2 


0.15 


0.023 


KCND3 


potassium voltage-gated channel, Shal-related subfamily, member 3 


-0.39 


< 0.001 


KCNE3 


potassium voltage-gated channel, Isk-related family, member 3 


0.23 


< 0.001 


KCNE4 


potassium voltage-gated channel, Isk-related family, member 4 


-0.38 


< 0.001 


KCNKl 


potassium channel, subfamily K, member 1 


0.25 


< 0.001 


KCNMAl 


potassium large conductance calcium-activated channel, subfamily M, alpha member 1 


-0.30 


< 0.001 


KCNN4 


potassium intermediate/small conductance calcium-activated channel, subfamily N, member 4 


0.23 


< 0.001 


MC0LN2 


mucolipin 2 


0.25 


< 0.001 


P2RX4 


purinergic receptor P2X, ligand-gated ion channel, 4 


-0.24 


< 0.001 


PKDl 


polycystic kidney disease 1 (autosomal dominant) 


-0.17 


0.012 


PKD2 


polycystic kidney disease 2 (autosomal dominant) 


-0.19 


0.004 


SCNIB 


sodium channel, voltage-gated, type 1, beta subunit 


-0.27 


< 0.001 


SCN7A 


sodium channel, voltage-gated, type VII, alpha subunit 


-0.41 


< 0.001 


SCNNIA 


sodium channel, non-voltage-gated 1 alpha subunit 


-0.18 


0.008 


TPCNl 


two pore segment channel 1 


-0.26 


< 0.001 


TPCN2 


two pore segment channel 2 


0.17 


0.013 


TRPCl 


transient receptor potential cation channel, subfamily C, member 1 


-0.25 


< 0.001 


TRPM4 


transient receptor potential cation channel, subfamily M, member 4 


-0.21 


0.002 


VDACl 


voltage-dependent anion channel 1 


0.21 


0.002 


VDAC2 


voltage-dependent anion channel 2 


0.22 


0.001 


VDAC3 


voltage-dependent anion channel 3 


0.26 


< 0.001 



is the Spearman's rank correlation coefficient. 
■^P-value Is calculated by Spearman's rank correlation test and adjusted by Benjamlnl & Hochberg correction. 



correlated with histological grade. There is a large over- 
lap between the three differentially expressed gene lists, 
which suggests a common ion channel-related biological 
mechanism underlying the different pathological pheno- 
types. The prognostic value of p53 mutation status 
has been well characterized in previous studies 
[15,29-31,35]. The frequency of p53 mutations is higher 
in ER negative breast cancer [25,35]. Mutant p53 and/or 
negative ER status are often associated with a high rate 
of proliferation, a high histological grade, and a poor 
prognosis [32,33,36,37]. Therefore, it is reasonable that 
several common ion channel genes are associated with 
p53 mutation status, ER status, and tumor histological 



grade, although the causal relationship between these 
three factors is still controversial. 

Gene expression-based molecular signatures have been 
proven as prognostically valuable in several human can- 
cers [38-42]. Gene signatures that work cooperatively 
with known clinicopathological factors may enhance 
prediction accuracy when identifying patients at higher 
risk for relapse and death. Our proposed molecular sig- 
nature that is composed of 30 ion channel genes (IC30) 
associated with tumor grade is a promising prognostic 
marker. IC30 was solely developed based on the discov- 
ery cohort and its prognostic power of IC30 was 
validated in seven independent validation cohorts. We 
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Figure 6 Kaplan-Meier curves for the patients from eight independent breast cancer cohorts. Panels from A to H denote the SIN, FRA, 
GER, NED, SWE, TWN, USAl, and USA2 cohorts, respectively. The expression of IC30 gene signature predicts poor survival in breast cancer. Red 
curves are for the IC30 positive patients while blue curves are for the IC30 negative patients. IC30 positive patients were defined as those having 
a IC30 risk score greater than the group median score. P-values were calculated by log-rank tests for the differences in survival between the IC30 
positive and negative groups. 



demonstrated that there were no significant multivariate 
interactions between IC30 status and other clinicopatho- 
logical covariates. When grouped by age, tumor size, 
tumor grade, or PR status, the expression of IC30 fur- 
ther stratified breast cancer patients with significant dif- 
ferences in survival A significantly increased risk of 
death was also observed in IC30 positive patients with 
positive lymph node status or positive ER status. How- 
ever, a significant difference between IC30 positive and 
negative groups among the patients with negative 
ER status was not detected, which may be due to the 

Table 5 Univariate Cox proportional hazards regression 
of survival by IC30 status for patients from 8 cohorts 



Cohort 


Hazard ratio 


95% Confidence interval 


P-value 


SIN 


1.98 


(1.15, 3.44) 


0.015 


FRA 


1.99 


(1.28, 3.10) 


0.002 


GER 


1.73 


(0.96, 3.14) 


0.069 


NED 


1.81 


(1.15, 2.86) 


0.010 


SWE 


4.33 


(1.76, 10.64) 


0.001 


TWN 


1.82 


(1.17, 2.85) 


0.008 


USAl 


3.11 


(2.05, 4.70) 


< 0.001 


USA2 


1.71 


(1.16, 2.51) 


0.006 



relatively smaller sample size in this category within the 
USAl cohort. In addition, we only detected a marginal 
significant association for patients with negative lymph 
node status in the USAl cohort, and a similar result was 
reproduced in the GER cohort. In fact, all patients from 
the GER cohort had negative lymph node status [34] 
and a marginally significant difference in the risk of 
death was observed between IC30 positive and negative 
patients in this cohort. Taken together with previous 
data, these results confirm that IC30 is not dependent 
on specific values of the respective covariates status, 
which enhances the identification of cancer patients at 
greater risk for death. 

Although IC30 gene signature was identified using the 
tumor grade information, both multivariate and univari- 
ate Cox regressions of survival reveal that IC30 is super- 
ior to tumor grade, which suggests that the prognostic 
information contained in IC30 is not limited to tumor 
grade. Tumor grade only explains a small proportion of 
variation (less than 25%) in expression of each IC30 
gene. The significant association between tumor grade 
and IC30 gene expression specifically implies that IC30 
genes are actively involved in tumor pathology in 
breast cancer. The combination of 30 genes contains 
more quantity of information than tumor grade itself. 
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Table 6 Multivariate Cox proportional hazard regression of survival for the patients from the USA1 cohort 



Covariat6 


HdZdrd rdtio 


95% Conficl6nc6 int6rval 


P-valu6 


IC30 + VS. - 


2.21 


(1.32, 3.70) 


0.002 


Age (per year) 


1.00 


(0.98, 1 .02) 


0.940 


Lymph node + vs. - 


2.07 


(1.35, 3.16) 


< 0.001 


Tumor size > T3 vs. < T3 


1.73 


(1.16, 2.57) 


0.007 


Grade 3 vs. 1,2 


0.67 


(0.41, 1.11) 


0.119 


ER + vs. - 


0.58 


(0.33, 1.01) 


0.055 


PR + vs. - 


0.76 


(0.45, 1.31) 


0.330 



Age < 50 



I N I m il 

^ ™l llllllllll M il 
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0 2 4 6 
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Figure 7 Kaplan-Meier curves for the patients grouped by clinicopathological factors. IC30 is independent from other clinicopathological 
factors in breast cancer. (A) Patients were stratified by age; (B) Patients were stratified by lymph node status; (C) Patients were stratified by tumor 
size; (D) Patients were stratified by tumor grade; (E) Patients were stratified by ER status; and (F) Patients were stratified by PR status. Red curves 
are for the IC30 positive patients while blue curves are for the IC30 negative patients. IC30 positive patients were defined as those having a IC30 
risk score greater than the group median score. P-values were calculated by log-rank tests for the differences in survival between the IC30 
positive and negative groups. 
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Table 7 Multivariate Cox proportional hazard regression of survival for the patients from the FRA cohort 


Covariate 


Hazard ratio 


95% Confidence interval 


P-value 


IC30 + VS. - 


2.55 


(1.21, 5.39) 


0.014 


Age (per year) 


1.00 


(0.98, 1 .02) 


0.687 


Grade 3 vs. 1,2 


0.79 


(0.43, 1 .45) 


0.444 


ER + vs. - 


0.91 


(0.33, 2.52) 


0.855 


PR + vs. - 


1.24 


(0.48, 3.18) 


0.655 


p53 mutant vs. wild-type 


1.22 


(0.70, 2.12) 


0.487 



Therefore, IC30 is based on tumor grade but not limited 
to tumor grade. 

The involvement of ion channels in human cancer has 
been intensively studied in the past years. However, 
there is no broad consensus on the role and interplay 
between ion channels and cancer. Generally, ion chan- 
nels are thought to "assist" cancer by tumor-related 
cellular behaviors such as proliferation, apoptosis, migra- 
tion, or angiogenesis [14,17,43,44]. However, it is diffi- 
cult to assign a detailed role for each ion channel in 
cancer pathology. In breast cancer, accumulating evi- 
dence indicates that channels play important roles in 
regulating tumor cell proliferation, cell cycle progression, 
and apoptosis [12,45,46]. Although a significant over- 
expression of channels has been correlated with 
human breast cancer cells [45,47-51], here we report a 
heterogeneous expression profiling in channel genes. 
Dysregulation in both directions was observed in 
channel genes within the IC30. Four channel genes 
in the IC30 were upregulated in high-grade tumors. For 
example, KCNN4, encoding an intermediate conduct- 
ance Kca channel (Kca3.1), is among the upregulated K"^ 
channel genes in the IC30 in high-grade tumors. The ex- 
pression pattern of Kca^.l was confirmed by a recently 
published study where Kca^.l mRNA and protein were 
more highly expressed in grade 3 tumors than in both 
grades 1 and 2 [52]. On the contrary, three K"^ channel 
genes within the IC30 were downregulated in tumors 
with higher grade, which includes KCNMAl encoding 
the BK channel alpha subunit. The negative correlation 
for KCNMAl between expression and tumor grade was 
in accord with four independent cohorts. However, in- 
creased expression of KCNMAl was found in metastatic 
breast cancer in the brain compared to metastatic breast 
cancers in other organs [53], which suggests a more 
complicated pathological role for the BK channel in 
tumor metastasis. 

Gene expression of CI" channels also demonstrated a 
heterogeneous pattern in breast cancer. We reported 
two up- and two downregulated CI" intracellular channel 
genes in high-grade tumors in the IC30. Among them, 
CLIC4 was found to be involved in skin cancer [54]; 
however, the exact role of CLIC4 is unclear. Besides CI" 
intracellular channels, the Ca"^ activated CI" channel 



CLCA2 was downregulated in breast cancer and is a 
candidate tumor suppressor gene [55]. We show here 
that the CLCA2 gene was upregulated in p53 mutant 
and/or ER negative breast tumors. In fact, the tumori- 
genicity of breast cancer was related with a loss of 
CLCA2 [56,57]. CLCA2 is a p53-inducible inhibitor of 
breast tumor proliferation [55]. However, the reason 
why CLCA2 expression is associated with p53 mutation 
status is beyond the scope of this study. 

is an essential regulator of the cell cycle and is indis- 
pensable for cell proliferation [14]. Increased expression of 
voltage-gated Ca^^ channels has been observed in colon 
cancer cells [58] and small cell lung cancers [59]. However, 
all 3 voltage-gated Ca^^ ion channel genes in the IC30 were 
downregulated in p53 mutant tumors and/or high-grade 
tumors. Apart from these voltage-gated Ca^"^ channels, the 
TRP family is also important to provide a Ca^^ influx path- 
way such that Ca^^ influx may occur through voltage-gated 
Ca^"^ channels [14]. Several studies have demonstrated that 
the expression of TRP channels is significantly upregulated 
in breast tumor tissue and breast cancer cell lines [60-63], 
which is related to malignant growth and cancer progres- 
sion [64]. However, a paradoxical result was observed in 
our study. PKDl, PKD2, and TRPCl, which are within the 
IC30 and belong to the TRP family, were downregulated in 
patients with p53 mutant tumors and/or of higher histo- 
logical grade. This expression pattern is consistent in sev- 
eral of the validation cohorts. 

Conflicting results were also seen in Na"^ channel genes. 
Decreased expression in high-grade tumors was found for 
the three Na^ channel genes in the IC30, which was con- 
firmed by the validation cohorts. However, increased ex- 
pression of voltage-gated Na^ channels has been reported 
in several cancer types, including breast, prostate, and lung 
cancer [65-68]. Na"^ channels were thought to enhance the 
invasiveness of cancer cells by increasing H"^ efflux [66] and 
by stimulating cysteine cathepsin activity [69]. The precise 
mechanism of Na^ channels in tumor development re- 
mains unclear [70]. The discrepancy between our results 
and previous observation may be due to the discrepancy 
between mRNA expression, protein expression, and chan- 
nel activity. Our study focused on mRNA abundance. 
However, protein expression and activity is not directly 
correlated to mRNA expression. Post-transcriptional 
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regulatory mechanisms predominantly control cellular 
mRNA to protein abundance ratios [71]. 

Voltage-dependent anion channels (VDACs) are a 
class of ion channel located in the outer mitochondrial 
membrane [72]. We observed a consistent and signifi- 
cant positive correlation between gene expression and 
tumor grade for the three VDAC genes in the IC30. A 
recently published study indicated that higher VDACl 
expression level predicts poor outcome in non-small cell 
lung cancers [73]. Here, we expanded this finding to 
breast cancer and the other 2 genes in VDAC family. 

Conclusions 

In summary, we investigated the gene expression profile 
of ion channels in breast cancer with respect to p53 mu- 
tation status, ER status, and histological grade. We show 
that there are numerous common ion channel genes, in- 
cluding ANOl, CACNAID, CACNA2D1, CACNA2D2, 
CLIC6, GLRB, KCND3, KCNE4, KCNMAl, KCNN4, 
MC0LN2, P2RX4, SCN7A, SCNNIA, and TPCNh that 
are differentially expressed with a change in p53 muta- 
tion status, ER status, and histological grade. The ex- 
pression pattern of some ion channels, including several 
potassium, calcium, and sodium channels, is contradict- 
ory to previously published results derived from breast 
cancer cell lines, animal models, and/or human patients. 
We also identified a molecular gene signature IC30, 
which represents a promising diagnostic and prognostic 
biomarker in breast cancer. Further investigation into 
the role of ion channels in tumor pathology could pro- 
vide new targets for therapy in multiple human cancers. 

Methods 

Ion channel genes 

The definition of human ion channel genes was obtained 
from lUPHAR-DB [74] and GeneCards [5,9]. In total, we 
collected 280 ion channel genes, including 5 Ca^ activated 
Cr channels, 6 CI' intracellular channels, 9 voltage- 
sensitive Cr channels, 1 mid-l-related CI' channel, 12 Kca 
channels, 48 Kv channels, 26 voltage-gated Ca"^ channels, 
14 NaV channels, 15 two-pore K"^ channels, 9 CatSper and 
two-pore channels, 16 inwardly rectifying K"^ channels, 4 
non-voltage-gated Na^ channels, 28 TRP channels, 10 cyc- 
lic nucleotide-regulated channels, 20 GABAa receptors, 5 
5-HT3 receptors, 5 glycine receptors, 18 ionotropic glu- 
tamate receptors, 16 nicotinic acetylcholine receptors, 7 
P2X receptors, 3 voltage-dependent anion channels, 1 
voltage-gated proton channel, 1 voltage-independent cat- 
ion channel, and 1 zinc activated ligand-gated ion channel 
(Additional file 1: Table SI). 

Gene expression data 

Eight independent microarray breast cancer datasets 
from Singapore (SIN) [15], France (ERA) [16], Germany 



(GER) [34], Netherlands (NED) [21], Sweden (SWE) 
[18], Taiwan (TWN) [19], and the United States (USAl 
and USA2) [20,27], were obtained for use in this study 
(Table 1). These datasets were chosen based on the large 
number of samples, the availability of clinical outcome 
data, and the diversity of tumor types. We assigned the 
SIN dataset as our discovery cohort and the other seven 
datasets as validation cohorts. 

Microarray data preprocessing 

The GC robust multichip average (GCRMA) algorithm 
[75] in Bioconductor was used to summarize the expres- 
sion level of each probe set for the microarray data from 
our discovery cohort (Affymetrix Human Genome U133 
set). Only the probe sets present (determined by func- 
tion "mas5calls" in the Bioconductor "affy" package) in 
at least one third of the samples were retained. We fur- 
ther limited our analysis to the probe sets with unique 
annotations and removed genes on chromosomes X and 
Y to avoid the potential confounding sex factor. 

Statistical analysis 

For the SIN and ERA cohorts, a two-tailed ^-test was 
used to identify the genes that were differentially 
expressed between p53 mutant and wild-type tumors. 
The genes with an adjusted P- value < 0.05 after 
Benjamini & Hochberg correction [22] and fold change > 
1.25 were deemed differentially expressed. The same 
methods and criteria were applied to identify the genes 
differentially expressed between ER positive and negative 
patients in SIN, ERA, USAl, and USA2 cohorts. 

The Spearmans rank correlation test was used to de- 
tect the relationship between ion channel gene expres- 
sion level and tumor histological grade. We calculated 
correlation coefficients and associated P-values using the 
R function "cor.test" with the "spearman" method. The 
genes with adjusted P- value < 0.05 after Benjamini & 
Hochberg correction were assigned as differentially 
expressed. We then tested the power of these tumor 
grade associated genes in predicting clinical outcome in 
breast cancers. Based on the relationship between gene 
expression and tumor grade, we can assign a Spearman s 
rank correlation coefficient to each gene as a weight 
(calculated solely from the discovery cohort). A risk 
score was then calculated for each patient using a linear 
combination of weighted gene expression as shown 
below: 

n 

i=l 

Here, s is the risk score of the patient; n is the number 
of differentially expressed genes; pi denotes the 
Spearmans rank correlation coefficient of gene /; Bi 



Ko et al. Molecular Cancer 201 3, 1 2:1 06 
http://www.nnolecular-cancer.conn/content/1 2/1/1 06 



Page 15 of 17 



denotes the expression level of gene /; and and T/ are 
the mean and standard deviation of the gene expression 
values for gene / across all samples, respectively. Patients 
were then divided into high-score (IC30 positive) and 
low-score (IC30 negative) groups with the median of the 
risk score as the threshold value. The median of the risk 
score was approximately equal to zero in each cohort 
(Additional file 2: Figure S5). A high score indicated a 
poor outcome. The "survival" library of the R was used 
to conduct survival analysis on the risk score. 

Hierarchical cluster analysis was conducted to gen- 
erate the gene expression heatmaps. The statistical 
significance of hierarchical cluster was evaluated by 
approximately unbiased P- value {AU), which is com- 
puted by multiscale bootstrap resampling. AU of 
cluster is a value between 0 and 1, which indicates 
how strong the cluster is supported by data. Higher 
AU means lower uncertainty of the hierarchical 
cluster. The "pvclust" library of the R was used to 
compute the AU values. 

Additional files 



Additional file 1: Table SI. Ion channel genes involved in this study. 
Table 52. Comparison in gene expression level between p53 mutant and 
wildtype tumors in validation cohorts. Table S3. Comparison in gene 
expression level between ER positive and negative tumors. Table S4. 
Correlation between gene expression and histological tumor grade. 
Table S5. Comparison in prognostic power between IC30 and 
clinicopathological factors for the USAl cohort. Hazard ratio was 
calculated separately for each variable by univariate Cox proportional 
hazard regression of survival. Table S6. Comparison in prognostic power 
between IC30 and clinicopathological factors for the FRA cohort. Hazard 
ratio was calculated separately for each variable by univariate Cox 
proportional hazard regression of survival. 

Additional file 2: Figure SI. Heatmaps of expression of the ion 
channel genes differentially expressed between ER positive and negative 
tumors. The differentially expressed genes were derived from the 
discovery cohort (SIN). Each row in the heatmaps was labelled with the 
corresponding gene symbol. The columns labelled with "-" denote ER 
positive tumors. Red represents relatively increased gene expression 
while blue represents down-regulation. Figure S2. Heatmaps of gene 
expression in FRA cohort. The listed genes are differentially expressed 
between ER positive and negative tumors in the discovery cohort (SIN). 
Each row in the heatmaps was labelled with the corresponding gene 
symbol. The columns labelled with "-" denote ER positive tumors. Red 
represents relatively increased gene expression while blue represents 
down-regulation. Figure S3. Heatmaps of gene expression in GER cohort. 
The listed genes are differentially expressed between ER positive and 
negative tumors in the discovery cohort (SIN). Each row in the heatmaps 
was labelled with the corresponding gene symbol. The columns labelled 
with "-" denote ER positive tumors. Red represents relatively increased 
gene expression while blue represents down-regulation. Figure S4. 
Heatmaps of gene expression in USAl cohort. The listed genes are 
differentially expressed between ER positive and negative tumors in the 
discovery cohort (SIN). Each row in the heatmaps was labelled with the 
corresponding gene symbol. The columns labelled with "-" denote ER 
positive tumors. Red represents relatively increased gene expression 
while blue represents down-regulation. Figure S5. Distribution of risk 
score. The red dash lines indicate the median of risk score. There is no 
significant deviation between zero and the median of risk score in each 
cohort (|z| < 0.2). 
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